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ABSTRACT 

We develop a time-dependent multi-group multidimensional relativistic radiative transfer code, 
which is required to numerically investigate radiation from relativistic fluids involved in, e.g., gamma- 
ray bursts and active galactic nuclei. The code is based on the spherical harmonic discrete ordinate 
method (SHDOM) that evaluates a source function including anisotropic scattering in spherical har¬ 
monics and implicitly solves the static radiative transfer equation with a ray tracing in discrete ordi¬ 
nates. We implement treatments of time dependence, multi-frequency bins, Lorentz transformation, 
and elastic Thomson and inelastic Compton scattering to the publicly available SHDOM code. Our 
code adopts a mixed frame approach; the source function is evaluated in the comoving frame whereas 
the radiative transfer equation is solved in the laboratory frame. This implementation is validated 
with various test problems and comparisons with results of a relativistic Monte Carlo code. These 
validations confirm that the code correctly calculates intensity and its evolution in the computational 
domain. The code enables us to obtain an Eddington tensor that relates first and third moments 
of intensity (energy density and radiation pressure) and is frequently used as a closure relation in 
radiation hydrodynamics calculations. 

Subject headings: relativistic processes — radiative transfer — shock waves 


1. INTRODUCTION 

Radiative transfer is an important piece of physics to 
describe how an astronomical object is observed. Also 
the radiation often affects dynamical behavior of the as¬ 
tronomical object. Thus it has been studied in many as¬ 
trophysical fields. However, the radiative transfer equa¬ 
tion is intrinsically a 6-dimensional Boltzmann equation 
that is computationally expensive. Therefore, the radia¬ 
tive transfer equation is frequently solved in a simplified 
and approximated form appropriate for an object of in¬ 
terest. 

Methods for the multidimensional radiative transfer 
and radiation hydrodynamics are developed in vari- 
ous fields, e.g., c osmol ogical structure formation (e.g., 
IHiev et all EiM [2?in9l for a review), star fo rmation 
(e.g.. iKrumholz et al.l 2?i?)7l : iTomid a et al.|| 2ni.‘I[l. a stel¬ 
lar and solar atmos phere (e.g., lAsnlund et al.l 1200(11 
Nordlund et al.ll2009ll. and a terrestrial atmosphere (e.g., 
Clough et al.ll2005HCollins et al.ll200^ . These methods 
are optimized for individual research fields and involve 
various sophisticated physics for individual phenomena, 
e.g., a chemical network for cosmological structure for¬ 
mation and star formation, fine-structure lines for a stel¬ 
lar and solar atmosphere, and molecular lines and scat¬ 
tering by dust for a terrestrial atmosphere. On the other 
hand, they ignore terms with higher orders than 0(y/c), 
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i.e., they assume that a fluid velocity is much slower than 
the light velocity. This is a commonly-used appropriate 
assumption to simplify the radiative transfer equation 
but inapplicable for the radiative transfer in a relativis¬ 
tic flow. 

An emission from a relativistic flow recently attracts 
researchers with the advent of a gamma-ray burst 
(CRB) 0 The CRB is a phenomenon emitting y-ray pho¬ 
tons from relativistic jets in a short period and one of the 
brightest objects in the Universe. The GRBs have been 
detected in t he distant Univers e with redshift a s high 
as z ~ 8.2 (jTanvir et al.l 120091 [Salrotgrra^e^jH J20W ; 
z ^ 9.4 for a photometric redshift. ICucchiara et al.ll2011 ) 
and are believed to probe the high-z Universe as well as 
quasars and galaxies. Furthermore, interestingly, it is ob- 
servationally exhibited that the y-ray emission has corre¬ 
lations of spectral peak en ergy versus isotropic radiation 
energy (lAmati et al.ll2002ll and spectral pea k energy ver¬ 
sus peak luminosity (|Yonetoku et all 120041 1. A correla¬ 
tion between X-ray luminosi ty at the break time and the 
break time is also suggested (jPainotti et al.l[2008ll . These 
correlations make researchers think of that the CRB can 
be a standardizable candle being detecta ble at higher 
redsh ift than Type la supernova (e.g., lAmati et al.l 

[2?iol . 

The mechanism of CRB prompt emission and the 
origin of correlations have been intensively studied. 
Observationally, many satellites and telescopes report 
large variations of the promp t emission, for ex ample, 
spectra with thermal (^c-g.- iR vde et all I20I0I1 ' non - 
thermal (e.g., lAbdo et al.l 1200911 iZhang fc Pe’edl2009H. 
and high- e nergy components (e.g.^ lAbdo et al.ll2009at 
Fan et all 1201311. polari zation (e.g., lYonetoku et al.l 
201lHUehara et al.ll2012ll . and duration (e.g., ultra-long 

® The emission from a relativistic flow is also interesting for, e.g., 
active galactic nuclei (AGN). 
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GRBs, iStratta et al.l 1201^ iLevan et al.l I20l4) . A cor¬ 
relation between optical and 7 -ray ligh t curves is also 
exhibited (e.g.. iVestrand et al.l 1200^ iWozniak et alJ 
120091: iGorbovskov et al.ll2012[l . Theoretically, the emis¬ 
sion mechan ism has been inv estigated by analytic stud¬ 
ies (e.g., IMeszaroi I2006L and references therein) 
or numerical studies with various assumptions: e.g., 
superposing blackbody radiation from a scattering 
photosphere ( Blinnikov et al.l 119991 : iMizuta et al.l 120111 : 
iNaeakura et al.l 120111 : iLazzati et ahl l2013[i n solving a 
radiative transfer e quation in a spherical steady flow 
(jBeloborodovI l20lHl . transferring photons in a steady 
flow with a relativ i stic Monte Carlo method (|G annio 


20001: IPe’ed 120081 iBeloborodovI 120111 : llto et ^ 


201 


Lnndman et al !r i201.H S. Shibata & N. Tominaga in 


prep.), and calc ulating a spherical relativistic rad iation 
hydrodynamics (|Tolstovll2005HTolstov et al.ll2013ll . 

In spite of plenty observations and theoretical stud¬ 
ies for many years, the mechanism of the GRB prompt 
emission is still under debate (e.g., lZhand[2014fl . This is 
mainly because most of studies is restricted to be qualita¬ 
tive and there are few quantitative studies taking account 
of a structure of relativistic jets. In order to investigate 
the GRB prompt emission quantitatively, a multidimen¬ 
sional relativistic radiation hydrodynamics calculation is 
essentially required. This is because the GRB is a rela¬ 
tivistic and multidimensional phenomenon and the radi¬ 
ation, of which energy can dominate the matter energy, 
closely couples with matter. Therefore, it is necessary to 
develop a radiative transfer code optimized for the GRBs 
which fully includes the terms higher than Oivjc). 

Much progress is recently made on the multidi¬ 
mensional radiation hydrodynamics calculation, for 
example, (I) special relativist ic 3-dimensional radi¬ 
ation magnetohydrodynam ics (iTakahashi et al.1 l2013t 
iTakahashi fc Ohsugal 1201311 . (2) genera l relativistic 3- 
dime nsional radiation hydrodynamics (|Sadowski et al.1 
nnn, m 3-dimensiona l radiation magnetohydrody- 
nami cs (l.liang et al.ll2012t iDavis et al.ll2012i l.liang et al.1 
[Ml . (4) relativistic Monte Garlo trans port coupled 
with hydrodynamics (jRoth fc KasenI 1^01511 . and (5) 3- 
di mensional special relati vistic Boltzmann hydrodynam¬ 
ics ([Nagakura et al.ll20H . However, the calculations (1) 
and (2) are based on the Ml closure method, which can 
treat anisotropic radiation field but not intersecting ra¬ 
diation from various sources. This is not suitable for 
the GRBs because the material, that is surrounding the 
relativistic jet, e.g., a cocoon, is hot and emits thermal 
photons to varying directions. The calculation (3) adopts 
the variable Eddington tensor (VET) method that can 
treat intersecting radiation from multiple sources and 
non-local radiation equilibrium, but ignores terms with 
orders higher than Oivjc). Recently. iJiang et al.1 (|20H 
implements a radiative transfer code involving time de¬ 
pendence and velocity dependent source terms but the 
method is still accurate up toOivjc). The calculation (4) 
couples a relativistic Monte Garlo transport method im¬ 
plicitly with hydrodynamics solvers and the calculation 
(5) is, in particular, adopted for the neutrino transport 
in a collapsing massive star. Although the calculations 
(4) and (5) would be applicable for the GRBs, the cal- 


^ The production site of photons is much deeper than the scat - 
tering photosphere in a relativistic flow (e.g., IShibata et al.ll2014ll . 


culation (4) can involve a noise of the Monte Garlo even 
with a reducing technique they developed and the cal¬ 
culation (5) is time-consuming and might be difficult to 
increase the number of mesh points because it involves 
an inversion of a huge matrix. There are also general 
relativistic radiative transfer codes to de scribe emission 
from surroundings of black holes (e.g., iDexter fc Agoll 
120091: [Shcherbakov fc Huandl2011h . The codes first de¬ 
rive photon geodesic trajectories in curved spacetime and 
integrate a relativistic radiative transfer equation along 
the geodesic paths. They also would be applicable for 
the GRBs. However, they waste time and computational 
resources because it is only required to calculate the tra¬ 
jectory in the flat spacetime in the GRBs. 

In this paper, we develop an implicit time-dependent 
multi-group multidimensional special relativistic radia¬ 
tive transfer (RRT) code with the mixed-frame approach. 
The RRT code can calculate an Eddington tensor with 
the intensity and thus could be the first step to a time- 
dependent special relativistic multidimensional multi¬ 
group radiation hydrodynamics code optimized for the 
GRBs. The RRT code is based on the spherical harmonic 
discrete ordinate method (SHDOM) and takes into ac¬ 
count ray tracing, time dependence, Lorentz transforma¬ 
tion, elastic Thomson and inelastic Gompton scattering. 
We present numerical test problems with the RRT code. 

We describe the method in Section and present test 
problems in Section (3] The results are summarized in 
Section [H 

2. METHOD 

The RRT code is based on Spheri cal Ha r monic 
Discrete Ordinate M ethod (SHDOM, lEvansI 119981 : 
iPincus fc Evansl[2009fF[ code which is a publicly-available 
static monochromatic radiation transfer code solving the 
following equation with the A iteration (Picard iteration) 
method. 

n • V/i,(s, n) = -x^(s, n)I„{s, n) -b n), (1) 

where s is a total path along a ray, n is direction of travel 
of the photon, and Xj/, and are an intensity, an 
extinction coefficient, and an emission coefficient at fre¬ 
quency v, respectively. Here, the extinction coefficient is 
the net absorption coefficient of a^+a^, where a^, and ai, 
are absorption and scattering coefficients, respectively. 
The SHDOM code is originally developed for radiative 
transfer in the terrestrial atmosphere and thus correctly 
treats anisotropic source terms including emission and 
scattering with spherical harmonics. Equation o is in¬ 
tegrated for each ray described with discrete ordinates {9, 
(jj) with short characteristic method (ray tracing). The 
SHDOM code transforms an intensity and a source func¬ 
tion between spherical harmonics and discrete ordinates 
at every time step. The scheme of the SHDOM code, 
including the transformation between spherical harmon¬ 
ics and discrete ordinates, the ray tracing, and so on, is 
com prehensively described and validated in lEvansI (|1998[ 1 
and iPincus fc Evaii^ (|20Q9ll . 

In order to apply the SHDOM code to a special rel¬ 
ativistic phenomenon, especially a GRB, the following 
processes and physics are necessary to be included: (I) 
time dependence because a time step is needed to be 

® http://nit.Colorado.edu/shdom.html 
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short enough to capture a relativistic fluid that moves 
as fast as the light speed, i.e., the light-crossing time is 
as long as the characteristic dynamical time, (2) Lorentz 
transformation between laboratory and comoving frames 
which results in relativistic beaming and a variation 
of photon frequencies, and (3) anisotropic and inelastic 
Compton scattering because the photon energy is com¬ 
parable with the electron rest-mass energy, thus the as¬ 
sumption of Thomson scattering is not valid, and the 
scattering dom inates the opacity in the relativistic jets 
of GRBs (e.g., lBeloborodovll20i3h . 

We set a computational domain to be translationally 
symmetric along the y axis but the ray of radiation is 
solved in three dimension and described with polar co¬ 
ordinates 6 and (j). The coordinates are set to have 
zenith direction of the z axis and azimuth angle mea¬ 
sured from the x axis to the y axis. The radiative transfer 
equation is solved with the mixed-frame approach (e.g., 
iMihalas fc Mihalaslll984 [Hubenv &: Burrov^ 120071) : the 

photon ray is traced in the laboratory frame and the 
source term is evaluated in the comoving frame. In 
addition to the above implementation, we update the 
ray tracing scheme to a cubic Bezier inte rpolant method 
(jde la Cruz Rodriguez fc Piskunovl[2013l l and include an 
accelerat ion scheme of the A iteration (Ng acceleration, 
iNd I1974D . but omit the adaptive treatment of mesh 
points and the parallelization with Message Passing In¬ 
terface (MPI) for simplicity in this paper. 


2.1. Time dependence 

A time-dependent radiative transfer equation is writ¬ 
ten as 

c 


dt 


+ n- s, n) 


= s, s, n) -b s, n), (2) 


where c is the light speed. A finite difference approxima¬ 
tion to the time derivative is written as 

dl^(t,s,n) I^(t + At,s,n) - I^{t,s,n) 

m “ At ■ 

The time-dependent radiative transfer equation is de¬ 
formed to the same shape as the static radiative transfer 
equation (Eq. [1} with modified absorption and emission 
coefficients as 


n • + At, s,n) = -Xvlu{t + At, s, n) -{-fj,,, (4) 

where 

= + At,s,n)-b ^ (5) 

~ / A \ ^ "I / ^\ 

= + At,s,n)-b———. (6) 

Since I^{t, s, n) and At are given, the equation (jU can be 
implicitly solved with the original SHDOM code. Such 
an implementation of time depende nce has been suc¬ 
cessfully adopted in previous studies (iBaron et al.l l2009t 
iHillier fc Dessartll2012t iJack et ^1201211 . ~ We note that 
the 4th order Runge-Kutta scheme is adopted to proceed 
the time step (Appendix lAT) . 


Variables at the next time step are evaluated in the lab¬ 
oratory frame with discrete ordinates with the ray trac¬ 
ing method. The variables are first Lorentz transformed 
with discrete ordinates to the comoving frame before the 
conversion from the discrete ordinates to the spherical 
harmonics. Then, the source terms are evaluated in the 
comoving frame with spherical harmonics. The source 
terms are first converted from spherical harmonics to dis¬ 
crete ordinates and then Lorentz transformed from the 
comoving frame to the laboratory frame. The obtained 
source term is adopted for the ray tracing in the labora¬ 
tory frame. 

The Lorentz transformation from the laboratory frame 
to the comoving frame moving with v, corresponding to 
the Lorentz fac tor of T, is defined with the following 
equations (e.g., IMihalas fc Mihal^ll984|l 


T. A 'n v 

vq = Tv[1 - 

V c 

V { V 

no = — 1 JT' — T— 
Vo { c 



n ■ V 


r-b 1 


( 7 ) 

( 8 ) 


where physical values in the comoving frame are denoted 
with a subscript 0. 

The Lorentz transformation is required for an inten¬ 
sity I, an emission coefficient rj, and an extinction coef¬ 
ficient Xi whose Lorentz invariants are I/v^, rj/v'^, and 
XV, respectively. The frequency and traveling direction 
of radiation are transformed with Equations © and ©a 
We prepare {0,(j)) mesh points in the laboratory and co¬ 
moving frames and the Lorentz invariants along the rays, 
which are transformed from a frame and remapped to the 
rays in the other frame, and then the intensity, emission 
coefficient, and absorption coefficient are obtained in the 
other frame. 


2.3. Thomson and Compton scattering 

The source function due to scattering is evaluated 
with spherical harmonics as the original S HDOM code 
does. Since the prescription is expatiated in lEvansI (jl993L 
119^ . we briefly describe the procedure. The source 
function (j)) is expanded in spherical harmonics space 
as 

S{fJ,,(j)) (9) 

Im 

where Yim{y,<p) sre orthonormal real-valued spherical 
harmonics functions, whereas the phase function of the 
scattering ^ is expanded in a Legendre series in the 
scattering angle as 


da 

dQ 


Nl 

^XiPh 


1=0 


( 10 ) 


where and Vi are the maximum order of Legen¬ 
dre polynomials, and Legendre polynomials, respectively. 
The source function is computed as 


Sim 


^Xi 

21 + 1 


Ilm + Ti 


m 1 


( 11 ) 


2.2. Lorentz transformation 


^ The subroutine is taken from 

http://cernlib.web.cern.ch/cernlib/version.html 
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Fig. 1. — Snapshots of mean intensity J for (a) t = 0.1 s and a = 0, (b) t = 0.3 s and a = 0, (c) t = 0.1 s and a = 10“^'^ and (d) 

t = 0.3 s and a = 10“^® cm“^. The medium is at rest (u = 0) and the x and z axes are divided by 512 mesh points. 


where w is the single scattering albedo, and Iim and Tim 
are intensity and thermal emission expanded in spherical 
harmonics, respectively. 

We adopt the phase functions of Thomson and Comp¬ 
ton scattering. The phase function of Thomson scatter¬ 
ing is 

^ = |(l + cos»e). (12) 


where tq is t he classical electron radius a nd 0 is the scat¬ 
tering angle (iRvbicki fc Lightmanlll985ll . Klein-Nishina 
scattering differential cross section is adopted for the 
Compton scattering. The equations are 




da 

dQ 


i + ;;^(i-cos0) 


I’d 


1^1 


( — -b -sin^0 
1^1 u 


(13) 

(14) 


where vi is the photon frequency after the scattering 
and t Ub is the rest mass of electron (IRvbicki fc LightmanI 

[HU). 

A change of the photon frequency and dependence of 
the scattering kernel on the photon frequency are essen¬ 
tial features of the Compton scattering. Therefore, we 


implement a multi-group treatment for test calculations 
of the Compton scattering. The differential cross section 
^ I (*)’(!) of an incident photon in i-th frequency bin with 
a range of -I- to j-th frequency bin with 


a range of 



,U) 


-b Av 


U) 

1 


is a frequency-dependent 


scattering kernel defined by Equation (fHl) . is 

expanded in a Legendre series in the scattering angle 
with frequency-dependent single scattering albedo. In 
the mixed-frame approach, the photon exchange between 
frequency bins takes place only in the comoving frame. 


3. TEST PROBLEMS 

lEvansI (1199811 and iPincus fc: EvansI (I2009f l have inten¬ 
sively tested the original SHDOM code and investigated 
its efficiency, accuracy, and scalability. Therefore, we 
focus on validation tests of time dependence, Lorentz 
transformation, and Thomson and Compton scattering 
in this paper. 


3.1. Searchlight beam test 

A searchlight b eam test has been performed in var¬ 
ious studies (e.g.. iRichling et al.l 120011: [Turner fc Stond 
120011 : iTakahashi et al.ll2013t iTakahashi fc Ohsugall2013ll . 
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Fig. 2.— Profile section at 2 : = 5.0 X 10® cm for the case of o = 0 
at t = 0.25 s (red) and an analytical result (gray). 

A narrow beam of light is introduced into the computa¬ 
tional domain at a certain angle and the beam crosses the 
domain. This test examines whether a code can solve a 
time-dependent radiative transfer equation and how ra¬ 
diation disperses along the path. 

We set a computational domain of 10^°cm square, of 
which X and z axes are divided by 512 mesh points. 
Numbers of angular mesh points are {Ng.N^) = (4,8). 
The origin is the bottom left corner of the domain. A 
beam of light is injected from the bottom boundary at 
1.5 X 10® cm < X < 4.5 x 10® cm along a ray with an 
angle of (6>, </>) = (O.ITtt, 0). Two cases of absorption co¬ 
efficients are tested; the domain is uniformly hlled with 
a medium with a = 0 or a = 10“^® cm“^. 

Figures[T](a)-[TKd) show snapshots of the mean intensity 
J in the domain at f = 0.1 s and 0.3 s. The beam crosses 
the domain properly with time at the injected angle. Fig¬ 
ure H] shows the prohle section at z = 5.0 x 10® cm and 
t = 0.25 s for the case of a = 0. The radiation disperses 
slightly laterally because the RRT code adopts the short 
characteristic method. 

Figures [3Ka)-[3Db) show J along x = tan (O.ITtt) z -|- 
3.0 X 10® cm. The Hgures demonstrate that a wave 
front proceeds with the light speed. However, the wave 
front is smeared. This is because the time dependence 
is taken into account with y and fj which exponentially 
reduce and/or enhance the intensity as a function of 
a path length s. The widths of the smooth prohle at 
the wave front are identical for the cases of a = 0 and 
a = 10“^® cm“^. The mean intensity of the beam is con¬ 
stant for the case of a = 0 and is reduced in accordance 
with exp (—as) for the case of a = 10“^® cm“^. The test 
conhrms that the time dependence and the radiation at¬ 
tenuation are correctly solved by the RRT code although 
the wave front is smeared. 

3.2. Two beam shadow test 

A s hadow test was proposed in iHaves fc Normaiil 
(|2003D . This test examines a reproduction of shadow 
behind an optically thick blob when plane-parallel radia¬ 
tion illuminates the blob. It is required to properly take 
into account at least up to the hrst moment of intensity 
in order to reproduce the shadow. Thus, for instance, a 
method with a diffusion approximation fails this test. 

On the other hand, a two beam test has be e n per ¬ 
formed in e.g., iDavis et’aTI (l201§) : iJiang et al.l (I2012D : 
iSadowski et al.l (|2013|) : iJiang et ail (|2014D . This test ex- 


Path length [c=1] 

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 



Fig. 3.— Time evolution of the mean intensity J along a ray 
with X = tan (O.ITtt) 2 + 3.0 X 10® cm at t = 0.05 s (red), 0.1 s 
(green), 0.15 s (blue), 0.2 s (magenta), 0.25 s (cyan), and 0.3 s 
(yellow) for the cases of (a) a = 0 and (b) a = 10“^® cm“^. The 
vertical dashed lines represent analytical expectations of the wave 
front proceeding with c at t = 0.05, 0.1, 0.15, 0.2, 0.25, and 0.3 s 
(gray dashed). The analytical result of radiation attenuation is 
also shown (black). 

amines whether two independent beams proceeding at 
different angles pass through without any interactions 
when they intersect. To describe this phenomenon, it is 
required to properly take into account at least up to the 
second moment of intensity. Thus, approximate schemes 
with a closure relation treating up to the hrst moment, 
e.g., the Ml closure method, fail this test. 

We solve a two beam with shado w test combining th e 
above tests , which is performed in iDavis et all (I2012D; 
[, iang et "ffil ()2ni2[ ): ISadowski et al.l ([20131 ): I.Tiang et all 
( 20141) . We set an optically thin (a = 0) computa¬ 
tional domain of 10^®cm square, where x and z axes 
are divided by 512 mesh points. Number of angular 
mesh points are {Ng,Nri)) = (4,8). The origin is the 
bottom left corner of the domain. An optically thick 
absorptive cylinder with a radius of 1.5 x 10®cm is lo¬ 
cated at (x, z) = (5.0 X 10® cm, 3.3 x 10® cm) per¬ 
pendicular to the xz plane. The optical depth of the 
cylinder is set to be r = 100 with the diameter, i.e., 
a = 3.3 X 10“® cm“^. Plane-parallel radiation is injected 
from the bottom boundary at angles of (0, (j)) = (O.ITtt, 0) 
and (O.ITtt, tt). 

Figures |4^-l4)a show snapshots of J in the domain at 
t = 0.2 and 1.0 s. The plane-parallel radiation proceeds 
properly with time and the wave speed is the speed of 
light. A shadow develops behind the cylinder along the 
directions of two beams and the two beams cross around 
(x, z) = (5.0 X 10® cm, T.4 x 10® cm) without any inter¬ 
action. The test conhrms that the calculation properly 
solves the radiation transfer equation for intensity. 

3.3. Radiative pulse test 

The RRT code is tested with the evolution of a ra¬ 
diative pulse initially having a Gaussian prohle in an 
optically thin medium and a scattering-dominated op- 
ticahy thick medium. The tests have been performed in 
ISadowski et al.l (I2013D and verify the treatments of time 
dependence and scattering. 
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Fig. 4.— Snapshots of the mean intensity J for the two beam and shadow test at (a) t = 0.2 s and (b) 1.0 s. The light is injected from 
the bottom boundary at the angles with (^, </>) = (O.ITtt, 0) and (O.ITtTjTt). An optically thick cylinder with a = 3.3 x 10“® cm~^ and a 
radius of 1.5 x 10^ cm is located at (x,z) = (5.0 x 10^ cm, 3.3 x lO^cm). The medium is at rest (u = 0) and the x and 2 : axes are divided 
by 512 mesh points. 
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Fig. 5.— (Top) Snapshots of the mean intensity J for the optically thin (y = 0) radiative pulse test at t = 0.1 s (left), 0.2 s (middle), 
and 0.3 s (right). (Bottom) Profile section of the mean intensity along the x axis (red), z axis (green), and x = z (blue) at t = 0.1 s (left), 
0.2 s (middle), and 0.3 s (right). The 1/r dependence is also shown in the bottom panels (magenta). 


3.3.1. Optically thin medium 

In the optically thin limit, the isotropic radiative pulse 
spreads with the speed of light and the mean intensity 
decreases inversely proportionally to the radius in the 
translational symmetry. 

We set an optically thin computational domain of 
10^'^cm square with a = 0. The origin is at the center 
of the domain. Each axis is divided by 128 mesh points. 
Number of angular mesh points are {Ng, N^) = (64,128). 
An initial radiative pulse is set according to the equation 


Io{x,n) = 100 exp 



(15) 


where r (= |a:|) is the distance from the origin and 
w = 9.0 X 10®cm. The intensity of the pulse is initially 


isotropic at each mesh point. 

The top panels of Figure [5] show snapshots of J at 
t = 0.1 s, 0.2 s, and 0.3 s. The expected size of the 
pulse is also shown by a green circle with a radius of 
ct. The bottom panels of Figure [S] show J along x axis, 
z axis, and a line with x = z aX t = 0.1 s, 0.2 s, and 
0.3 s. These profiles of J are identical and the peak of 
J decreases with the expansion of the pulse according 
to J oc 1/r as expected. The result demonstrates that 
the pulse isotropically propagates with the speed of light 
and that the geometrical dilution of radiation is correctly 
followed with the RRT code. 

3.3.2. Optically thick medium 

In the scattering-dominated optically thick medium, 
the radiative pulse diffuses out. A one-dimensional dif- 




























7 



-0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 

X[10^°cm] 


Fig. 6.— Profile section of the mean intensity J along the x axis 
(crosses) for the optically thick (cr = 3 x 10~® cm“^) radiative 
pulse test at t = 0 (red), 0.5 s (green), 1.0 s (blue), and 3.0 s 
(magenta). The solid lines show the exact solution (Eg. HOD . 

fusion equation ^ can be solved analytically 

and the solution is written as 

( 16 ) 

where x = ^-^d Jo{x) is the initial profile of the 

radiative pulse. 

We set the optically thick computational domain of 
10^°cm square with r = 300 along a side, i.e., ct = 3 x 
10“® cm“^. The origin is set at the center of the domain. 
Each axis is divided by 128 mesh points. Number of 
angular mesh points are {Ng,N^) = (32,64). The each 
mesh point is optically thick, r = 2.3, so that most of 
the photons are scattered at least once in the mesh point. 
Here, we adopt a spherical scattering kernel. A radiative 
pulse is initially set according to 


Io{x,n) = 100 exp 



(17) 


where w = 9.0 x 10®cm. The intensity of the pulse is 
initially isotropic at each mesh point. 

Figure [6] shows the time evolution of J along the x axis. 
The radiation gradually diffuses in the numerical simula¬ 
tion and the distribution of J well reproduces that of the 
analytic solution at every time step until the radiation 
reaches the boundaries of the computational domain. 


3.4. Relativistic beaming test 

We validate the ability of the RRT code to treat the 
relativistic beaming as a result of the Lorentz trans¬ 
formation. This is a characteristic feature of a special 
relativistic radiative transfer calculation. Although the 
propagation of beam of light along a curved trajectory 
is demonstrated i n general relativistic r adiative transfer 
calculations (e.g., iSadowski et al.l[^013f) . this is the first 
attempt to confirm the relativistic beaming in special 
relativistic calculations. 

We set a computational domain of 10^° cm square filled 
with the optically thin medium with a = 0 and consider 
a cylindrical light source with a radius of i? = 5 x 10®cm 
at the center of a computational domain. The center of 
a computational domain is set to be the origin. In the 
cylinder, particles move with v = (vx,Vz) with respect 
to the laboratory frame and emit photons isotropically in 


the comoving frame of each particle, whereas they do not 
emit any photons and not interact with photons outside 
the cylinder. Each axis is divided by 128 mesh points and 
number of angular mesh points are {Ng,N^) = (32,64). 

Figures [7^-[7f show snapshots of the normalized mean 
intensity J/Jmax in the Laboratory frame at t = 0.2 s 
for the cases with {vx,Vz) = (0,0), (0.1c, 0), (0.5c, 0), 
(0.9c, 0), (0.99c, 0), and (0.995c, 0), respectively, where 
Jmax is the maximum mean intensity in the computa¬ 
tional domain. The corresponding Lorentz factors are 
r = 1, 1.005, 1.15, 2.29, 7.09, and 10.0. It is clearly 
shown that the mean intensity is high along the direc¬ 
tion of V. An analytic expression of the beaming effect 
(asymptotically 9 ^ jr- l^'^bicki fc LightmanlflOSSll is also 
shown with green lines circumscribing the light source. 
The biased distribution of J in the calculation is consis¬ 
tent with the analytic expression. The test shows that 
the beaming effect is correctly taken into account in the 
RRT code. 

The speed of 0.995c is close to the maximum speed 
that can be correctly solved with Ng = 32 because the 
half-angle of the concentration of radiation for the case 
of u = 0.9952c is comparable to the interval between 
angular mesh points. Although the test is limited due 
to computational resources, the Lorentz factor as high 
as the one encountered in GRB models can be resolved 
by the RRT code, if the larger number of angular mesh 
points is adopted, because the mapping subroutine for 
the Lorentz transformation is valid even for high Lorentz 
factor. Figure [5] shows the angular distribution of the 
intensity I in the laboratory frame for the cases with 
ivx,c-Vz) = (0, 5 X lO-^c) (T = 10), (0, 5 x IQ-^c) (T = 
100), (0,6 X lO-^c) (T = 300), and (0,5 x W^c) (T = 
1000), which is transformed from an isotropic intensity Iq 
in the comoving frame. The deviation from the analytical 
solution J(^) = (r(l — fiVz/c))~^ Iq is less than 10“®/. 
Here, we adopt the large number of angular mesh points 
of {Ng,N^) = (4096,8192). We note that the mapping 
subroutine correctly transforms an isotropic radiation for 
any Lorentz factor even with the small number of angular 
mesh points. 

3.5. Comparison with Monte Carlo method 

The original SHDOM code and the Monte Carlo 
method have been carefully compared an d their draw¬ 
bac k and advantage have b een presented in lEvansI (119981) 
and iPincus fc EvaS (|2009ll . Therefore, we test the RRT 
code only on the treatments of the Thomson and Comp¬ 
ton scattering and Lorentz transformation, and the im¬ 
plementation of multiple frequency groups for the Comp¬ 
ton scattering tests, by comparing solutions of a shadow 
test of the RRT code with those of a relativistic Monte 
Carlo (RMC) code (S. Shibata & N. Tominaga in prep.. 
Appendix B for a test of the RMC code). In this subsec¬ 
tion, we omit the time dependence and adopt the static 
version of the RRT code because the RMC code currently 
does not follow the time evolution of the radiation. 

Here, we adopt an ideal test problem. A 
scattering-dominated optically thick cylinder with cr = 

tTo exp (—^), where cto = 2 x 10“^ cm“^ and re = 2.5 x 

10® cm, is put at the center of the optically thin computa¬ 
tional domain of 10^® cm square with a = 2x 10“^® cm“^. 
Scattering particles rest or flow with v in the scatter- 
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Fig. 7.— Snapshots of the mean intensity normalized with the maximum mean intensity J/ Jm-c 
the cases of (a) (vx,Vz) = (0,0), (b) (vx,Vz) = (0.1c,0) (c) {vx,Vz) = (0.5c,0), (d) (va:,Vz) = 
{yx,Vz) = (0.995c, 0). Asymptotic analytic expressions of the beaming effect are also shown (green 



Fig. 8.— Distribution of the normalized intensity I/Io as a func¬ 
tion of the zenith angle for the numerical solution (red) and the 
analytic solution (green). 

ing cylinder^ Plane-parallel radiation is injected from 
the left and bottom boundaries at a single angle of 
{0,(j)) = (0.247r,0). The incident angle of radiation is 
set to the same in the RRT and RMC calculations. 

3 . 5 . 1 . Thomson scattering 

In this subsection, we adopt a kernel of Thomson scat¬ 
tering. The computational domain is divided by 128 x 128 
mesh points and numbers of angular mesh points are 
{Ne,N^) = (160,320). 

Figures [i(a)-ini;b), and [ni[c)-|nild) show normalized z- 
oriented fluxes Fz{x, z)/Fz^ma.^ obtained by the RRT and 
RMC calculations in the cases with (vx.Vz) = (0,0) and 

We assume, for example, an ideal atom that do not interact 
with the photons if they are neutral and that the atoms rest or 
flow with V in the computational domain and its ionization fraction 
changes to reproduce the distribution of a. 


for the relativistic beaming test for 
(0.9c,0), (e) {vx,Vz) = (0.99c,0), and (f) 
lines). 

(0,0.9c), respectively, where is the maximum z- 

oriented flux in the computational domain. Figures IHKa)- 
[5Kb) demonstrate that the z-oriented fluxes are small be¬ 
low the cylinder due to the scattered photons propagat¬ 
ing to the — z direction, whereas the z-oriented fluxes are 
high at the top left side of the cylinder. Since the number 
density of scattered photons decreases with the distance 
from the scattering cylinder as 1/r, the z-oriented fluxes 
are prominently modified around the cylinder. Also the 
shadow behind the cylinder consistently forms in both 
calculations. Figures inKc)-|5Kd) shows that the scattered 
photons are concentrated to the +z direction due to the 
relativistic beaming in both calculations. These demon¬ 
strate that the RRT code well solves the anisotropic scat¬ 
tering and the Lorentz transformation. 

3 . 5 . 2 . Compton scattering 

In this subsection, we adopt frequency-dependent ker¬ 
nels of Compton scattering and the multi-group treat¬ 
ment. The computational domain is divided by 64 x 64 
mesh points. We test the RRT code for 2 cases with 
different velocity of scattering particles with {vx,Vz) = 
(0,0) and (0,0.9c). Numbers of angular mesh points are 
{Ng,Nff,) = (128,256). The frequency range is equally 
divided by 10 bins in both cases but the maximum fre¬ 
quency is hv = l.lrUeC^ for the cases with {vx,Vz) = 
(0,0) and 4.0meC^ for the case with {vx,Vz) = (0,0.9c). 
Monochromatic light with hiy = l.OSweC^ for the cases 
with (vx, Vz) = (0,0) and hv = l.OrUeC^ for the case with 
{vx,Vz) = (0,0.9c) is injected. 

Figures IHaj-IIOKf) and [iira)p:f) show 
frequency-dependent normalized z-oriented fluxes 
Fu,z{x,z)/F^^z,max in the cases with {vx,Vz) = (0,0) and 
(0,0.9c), respectively, where F^^z,m&x is the maximum 
z-oriented flux in a frequency bin in the computational 
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Fig. 9.— Snapshots of 2 ;-oriented fluxes normalized with the maximum flux Fz{x, z)/Fz^max in the medium with {vx,Vz) = (0,0) [(a) 
and (b)] and {vx,Vz) = (0, 0.9c) [(c) and (d)]. The panels show the results of the RRT code [(a) and (c)] and the RMC code [(b) and (d)]. 
Here, the Thomson scattering kernel is adopted and the plane-parallel radiation is injected from the left and bottom boundaries at a single 
angle of {6, <f>) = (0.247r, 0). 


domain. Each panel shows the fluxes in different energy 
bins. There was no light in these energy bins before 
scattering. This is the reason why no shadow appears 
behind the cylinder. The negative fluxes due to the 
scattered photons produce the black region below the 
cylinder in Figure [TOT el and [TOT f). Although the shape 
of F^^z{x,z) is smeared in the RRT code, especially in 
the case of {vx,Vz) = (0,0.9c), because the frequency of 
photons is converted to the central frequency of each bin 
at every time step, it provides the similar snapshots of 
the z-oriented flux in the RMC code. These display that 
the RRT code correctly solves the Compton scattering 
and the Lorentz transformation, and well treats the 
multi-group radiation transfer equation. 

4. SUMMARY 

We develop a time-dependent multi-group multidimen¬ 
sional relativistic radiative transfer code by implement¬ 
ing the treatments of time dependence, multi-frequency 
bins, Lorentz transformation, and elastic Thomson and 
inelastic Compton scattering to the publicly available 
SHDOM code. The SHDOM code evaluates a source 


function in spherical harmonics and solves a static ra¬ 
diative transfer equation with a ray tracing in discrete 
ordinates. The RRT code is validated by the various 
tests and the comparison with the RMC calculations. 
The searchlight beam, two beam with shadow, radia¬ 
tive pulse, and relativistic beaming tests are successfully 
passed by the RRT code and confirm that the RRT code 
correctly handles the time dependence and the Lorentz 
transformation. The results of the RRT code are consis¬ 
tent with those of the RMC code and the comparisons 
verify the implementation of elastic Thomson and inelas¬ 
tic Compton scattering and multi-group treatment in the 
RRT code. The RMC code, in turn, is validated against 
the ECS tools (Appendix [Q. 

The RRT code enables us to obtain the evolution of 
intensity and thus to self-consistently derive an Edding¬ 
ton tensors without approximations like the flux lim¬ 
ited diffusion or the Ml closure methods. We empha¬ 
size that the radiation tends to be more anisotropic in 
the relativistic fluid because of the Lorentz transforma¬ 
tion and Compton scattering for the y-ray photons, and 
thus that the angular distribution of radiation should 
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Fig. 10.— Snapshots of frequency-dependent 2 -oriented fluxes normalized with the maximum flux z)lFi,^z,ma^ in the cases with 

{vxiVz) = (0,0). The panels represent fluxes at the energy bins with [O.SSmeC^, 0.99meC^] [(a) and (b)], [0.66meC^, 0.77meC^] [(c) and 
(d)], and [0.44meC^, 0.55meC^] [(e) and (f)]. The panels show the results of the RRT code [(a), (c), and (e)] and the RMC code [(b), (d), 
and (f)]. Here, the frequency-dependent Compton scattering kernels are adopted. 


be properly taken into account. Combining the Edding¬ 
ton te nsors with relativ istic hydrodynamics calculations 
(e.g., iTominag^ 12009^ . a relativistic radiation hydro¬ 
dynamics will be realized with the variable Eddington 
tensor method. Furthermore, the RRT code implicitly 
solves the radiative transfer equation and thus can fol¬ 
low the radiative transfer in a non-relativistic fluid like 
a supernova without adopting unnecessarily short time 
steps. Such a method will be useful to clarify the con¬ 
nection between GRBs and supernovae. 

It is currently difficult to increase the numbers of fre¬ 
quency bins, and angular and spatial mesh points due to 
an available memory resource. The difficulties are solved 
if the adaptive treatment of mesh points and the paral¬ 
lelization with distributed memory are implemented be¬ 


cause the ray tracing in the laboratory frame and the 
evaluation of source function in the comoving frame are 
independent of each ray and each mesh point, respec¬ 
tively. However, we note the low resolution of the fre¬ 
quency is an intrinsic drawback of the multi-group treat¬ 
ment compared to the Monte Carlo method, in which a 
frequency of each photon changes continuously. Thus, 
the Monte Carlo method is superior to the RRT code 
for spectral synthesis calculations of lines and fine spec¬ 
tral features. However, the method intrinsically involves 
a noise and it is expensive to reduce the noise because 
the reduction is realized only proportional to the square 
root of the number of photon packetsE3 And, the large 

Several techniques are suggested to reduce the noise (e.g., 
ISteinacker et aLll2ni.^ : IRoth Kr, Kasenll2ni,^ l. 
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Fig. 11.— Same as Figure [To] but for the cases with {vx,Vz) = (0,0.9c). The panels represent fluxes at the energy bins with 
[2.0meC^, 2.4meC^] [(a) and (b)], [l.2meC^, l.GrrieC^] [(c) and (d)], and [0.4meC^, O.SmeC^] [(e) and (f)]. 


number of photon packets is necessary to follow the time 
dependence because the photon packets are emitted at 
each time step. Furthermore, the radiation contributes 
to the hydrodynamics with an integration of radiation 
over all frequency and the dynamical effects of the ra¬ 
diation in GRBs are not dominated by narrow spectral 
lines. Thus, we propose a post-processing Monte Carlo 
calculation for spectral synthesis after a time-dependent 
matter-coupled relativistic radiative transfer calculation. 
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APPENDIX 
TIME INTEGRATION 

We treat the time dependence with the modified absorption and emission coefficients as shown in Section 12.11 and 
implicitly obtain the self-consistent time-dependent intensity at t -|- At. Although the time derivative is differentiated 
in the first order, we adopt the 4th order Runge-Kutta scheme to proceed the time step from t to t At by dividing 
the time interval to 4 steps, in order to increase the accuracy by adopting smaller time intervals. 
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Fig. 12.— Fluence as a function of the depth calculated with the RMC code (thick) and the EGS tool (thin). The color represents the 
rings with r < 1 x 10^ cm (red), 1 x 10^ cm < r < 5 x 10^ cm (blue), 5 x 10^ cm < r < 1 x 10^ cm (black), and 1 x 10^ cm < r < 2 x 10^ cm 
(green). 


I,{t + At, s, n) = s, n) + + A/'^) + AJ^^) ^ ^j(i)^ 

where 

AI^"> = y,,{At/6,I^(t,s,n)) - Ii,{t,s,n) (A2) 

/«=4(l,s,n)+3A/« (A3) 

A/(2)=^,(AV3,/«)-/W (A4) 

(A5) 

A43) ^ ^^(Al/3,42)) - /(2) (A6) 

/(3)=/,(i,s,n)+3A/(3) (A7) 

A/(4)=./,(AV6,/(3))-/(3). (A8) 


Here, J^v{At, 1^) is the solution of Equation dU with 1,^ and At, i.e., J^i,(At, I^{t, s, n)) = Iy{t + At, s, n). 

COMPARISON WITH ELECTRON GAMMA SHOWER (EGS) SOETWARE 

We compare the result of the RMC code with that of National Research Council’s electron gamma shower (EGS) 
software tooQ to confirm its validity, especially the treatment of the scattering. EGS software is a publicly available 
code that enables sophisticated treatment of photon, electron, and positron transfer in complicated medium and 
provides graphic tools to show numerical results. Here, we adopt Rurznrc package only with the Compton scattering 
by electrons at rest in the EGS software. 

We set a cylindrical computational domain with a radius of 2 x 10® cm and a depth of 2 x 10® cm. The cylinder is 
filled with H atoms with a number density of 8.37 x 10“® g cm“®. The incident photons with 50 keV are vertically 
injected from the top boundary with r < 1 x 10® cm, where r is the distance from the center of the top boundary. 

Figure [T^ shows a comparison with the fluence as a function of the depth in rings with r < 1 x 10^ cm, 1 x 10"* cm < 
r < 5 X 10"^ cm, 5 x 10"* cm < r < 1 x 10® cm, and 1 x 10® cm < r < 2 x 10® cm, normalized by the maximum fluence. 
The fluence is defined as a summation of 1/| cos0| for all photons per unit area at each boundary of slabs, where 0 is 
the angle the photon makes with respect to the normal to the plane. We adopt 6 x 10^ photons for both of the RMC 
calculation and the calculation by the EGS tool. The RMC codes give consistent results with that obtained by the 

http://www.nrc-cnrc.gc.ca/eng/solutlons/advisory/egsnrc_index.html 
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EGS software. 
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